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1.  Introduction 


The  behavior  of  polymers  under  extreme  conditions  (high  pressure  and  temperature)  is  of  interest 
for  both  civilian  and  military  applications,  such  as  polymer-bonded  explosives,  coatings, 
adhesives,  and  lightweight  armor.  Insight  into  the  response  and  energy  dissipation  mechanisms 
of  such  materials  to  ballistic  impact  of  projectiles  is  necessary  in  order  to  design  and  develop 
novel  materials  to  improve  the  protection  of  the  Future  Force.  The  material  properties  and 
response  at  extreme  conditions  can  be  determined  through  shock  experiments,  which  are  often 
difficult  to  measure  experimentally  because  of  difficulties  in  traversing  a  large  range  of  pressures 
(up  to  hundreds  of  gigapascals)  and  temperatures  (thousands  of  kelvin)  with  available 
instrumentation.  In  addition,  interesting  behavior,  such  as  observed  behind  a  shock  front,  occurs 
at  extremely  short  time  and  length  scales  (nanoscale),  which  poses  problems  in  characterizing  the 
material  using  current  experimental  capabilities.  To  further  understand  the  shocked  systems, 
simulation  methods  such  as  molecular  dynamics  (MD)  and  quantum  mechanics  can  be  used  to 
provide  insight  into  atomic-level  phenomena. 

The  shock  response  of  a  material  can  be  described  by  its  shock  Hugoniot,  which  is  the  locus  of 
thermodynamical  states  accessible  by  shock  loading  from  a  given  thermodynamic  initial 
condition.  This  calculation  can  be  performed  through  a  variety  of  methods,  such  as  quantum 
mechanical  theories  (wave  function-based  methods  and  density  functional  theory)  and  classical 
atomistic  simulation  methods,  such  as  MD  and  reactive  Monte  Carlo.  Insight  into  atomic-level 
phenomenon  of  these  materials  at  extreme  conditions  is  most  readily  extracted  using  the  MD 
simulation  method,  where  regimes  not  accessible  by  experimental  techniques  can  be  accessed. 
MD  computation  of  the  Hugoniot  state  can  be  obtained  through  several  different  types  of 
simulations.  One  can  involve  the  calculation  of  properties  behind  the  shock  discontinuity  in  a 
shock  wave  simulation  (7).  Another  method  introduced  by  Erpenbeck  (2)  involves  generating  an 
equation  of  state  for  the  subsequent  evaluation  of  Hugoniot  conservation  relations  (5).  An 
equilibrium  uniaxial  Hugoniostat  method  can  also  be  employed  that  uses  equations  of  motion 
which  restrain  the  system  during  the  simulation  so  that  the  time  average  properties  amount  to  the 
Hugoniot  curve  ( 4 ). 

One  of  the  main  limitations  of  MD  is  that  the  accuracy  of  the  results  is  highly  dependent  on  the 
description  of  the  forces  acting  on  the  particles.  An  erroneous  description  can  lead  to  inaccurate 
results,  thus  it  is  important  to  determine  the  quality  and  applicability  of  the  force  field  that 
describes  the  interactions  between  molecules.  Therefore,  in  order  to  properly  depict  the  material 
properties  of  materials  subjected  to  the  extreme  temperatures  and  pressures  corresponding  to 
shock  conditions,  force  fields  must  also  be  accurate  for  both  the  shocked  and  unshocked  states. 

In  addition,  during  shock  experiments,  bonds  between  atoms  in  the  material  may  break  and  the 
formation  of  new  reaction  products  can  occur.  This  behavior  will  require  the  use  of  a  reactive 
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potential  (i.e.,  one  that  describes  making  and  breaking  of  chemical  bonds)  whose  parameterization 
is  currently  an  area  of  interest  (5).  Most  reactive  potentials  available  for  polymers  are  highly 
idealized  representations,  although  ReaxFF  is  emerging  as  one  that  can  properly  describe 
chemistry.  In  lieu  of  accurate  reactive  potentials,  numerous  studies  rely  on  nonreactive  potentials 
that  have  been  extensively  parameterized  and  tested  against  systems  at  ambient  conditions. 

The  purpose  of  this  study  is  to  assess  the  quality  of  a  popular  nonreactive  interaction  potential  in 
describing  the  shock  properties  of  amorphous  atactic  poly[methyl  methacrylate]  (PMMA) 
through  MD  simulation.  In  this  study,  we  use  the  polymer-consistent  force  field  (PCFF)  with 
condensed-phase  optimized  molecular  potentials  for  atomic  simulation  studies  (COMPASS) 
charges.  We  will  investigate  the  effect  of  varying  system  size,  equilibration  time,  and  chain 
length  on  the  calculated  Hugoniot  curve.  We  will  also  consider  relatively  low  pressures 
(moderate  compressed  systems)  where  bond  breaking  does  not  occur,  which  justifies  the  use  of  a 
nonreactive  force  field. 


2.  Methodology:  Hugoniot  Calculation 


Although  several  methods  exist  to  calculate  the  Hugoniot  curve,  we  evaluated  the  Hugoniot 
points  through  a  procedure  developed  by  Erpenbeck  (2),  which  involves  performing  several 
constant  particle,  volume,  and  temperature  (NVT)  simulations  at  multiple  temperatures  for 
several  compressed  structures.  The  Hugoniot  curve  consists  of  the  set  of  pressure-volume- 
temperature  points  for  which  the  Hugoniot  expression 

Hg  =  E  -  Eo  +  lMP  +  Po)(  V-  Vo)  (1) 

is  0.  In  this  equation,  E  is  the  specific  internal  energy  per  unit  mass  (sum  of  the  kinetic  and 
potential  energy),  P  is  the  pressure,  and  V  =  1/p  is  the  specific  volume  (p  is  the  density).  The 
term  specific  refers  to  the  quantity  per  unit  mass,  while  the  subscript  “o”  refers  to  the  quantity  in 
the  initial  unshocked  state. 

The  Hugoniot  points  are  calculated  through  several  MD  simulations,  where  a  series  of  NVT 
simulations  are  performed  over  a  range  of  temperatures  at  a  fixed  specific  volume  after  the 
system  is  annealed  and  relaxed.  To  obtain  systems  at  different  specific  volumes,  the  polymer  is 
compressed  isotropically  at  a  specified  pressure  and  allowed  to  equilibrate  under  constant 
temperature  and  pressure  (NPT)  MD  simulations  (T  =  298  K).  Using  the  equilibrated  structure 
obtained  at  the  desired  pressure,  we  performed  a  series  of  sequential  NVT-MD  simulations 
where  the  temperature  was  incrementally  increased  after  a  prescribed  simulation  time.  After  the 
temperature  scan,  the  Hugoniot  function  (equation  1)  is  evaluated  at  each  temperature;  this 
produces  a  series  of  equation-of-state  points  as  a  function  of  temperature  for  each  specific 
volume,  which  can  be  interpolated  through  a  linear  line  to  locate  the  temperature  (the  Hugoniot 
temperature)  at  which  the  expression  is  0.  Finally,  to  obtain  the  Hugoniot  pressure,  the 
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pressure  and  temperature  data  are  fitted  with  a  linear  equation  and  evaluated  at  the  Hugoniot 
temperature.  Since  the  linear  fit  was  only  performed  on  the  points  which  bracket  the  Hugoniot 
temperature  and  pressure,  the  temperature  scan  ceased  once  the  Hugoniot  point  was  bracketed. 
Further  NVT-MD  simulations  can  be  performed  to  ensure  that  these  two  points  are  at 
equilibration.  For  a  schematic,  see  figure  1 . 


Figure  1.  Schematic  of  Erpenbeck  method  used  to  calculate  Hugoniot  curve. 


Experimentalists  tend  to  report  shock  Hugoniot  data  in  terms  of  the  Ust  and  Upt ,  which  are  the 
shock  and  particle  (or  mass)  velocities,  respectively.  To  obtain  these  values  from  the  P  and  V 
data  points,  the  following  equations  can  be  used: 


U 


st 


WjzEia 

]  \-V!V0  ’ 


(2) 


and 

u„=pAP-PM-r/r,). 


(3) 


On  the  same  note,  the  experimental  data  can  be  converted  into  P  and  V  data  by  rearranging 
equations  2  and  3: 


and 


V  = 


U,„ 


p  _  p  _|_  w  U p{ 


V, 


(4) 

(5) 
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3.  Computational  Methods 


MD  simulations  were  performed  using  the  open-source  classical  molecular  dynamics  code  large- 
scale  atomic/molecular  massively  parallel  simulator  (LAMMPS)  ( 6 ,  7).  We  have  interfaced 
amorphous  polymer  builders  with  LAMMPS  using  the  commercial  visualization  package 
materials  processes  and  simulations  (MAPS)  (5).  MAPS  was  used  to  create  a  periodic  cell  of 
PMMA  and  assign  the  PCFF  force  field,  where  the  COMPASS  charges  were  obtained  using 
Materials  Studio  (9).  To  produce  the  three-dimensional  periodic  cell,  we  used  a  well-established 
method  to  build  the  so-called  amorphous  cells.  This  method  uses  a  Monte-Carlo  technique  to 
build  an  amorphous  structure  of  the  polymer  followed  by  the  energy  minimization  of  the  bulk 
structure.  The  scheme  builds  polymeric  chains  by  using  rotational  isomeric  state  theory  and 
taking  into  account  nonbonded  interactions  with  the  already  constructed  neighboring  chains, 
including  their  periodic  images  (JO). 

PCFF  is  a  class  2  force  field  where  the  nonbonded  interactions  are  composed  of  a  9-6  Lennard 
Jones  potential  (van  der  Waals)  and  a  Coulombic  pairwise  (electrostatic)  interaction.  In  our 
simulation,  we  used  a  cutoff  of  12.0  A  for  both  interactions.  We  also  included  an  accurate 
description  of  the  long-range  electrostatic  interaction  using  the  particle-particle,  particle-mesh 
method  with  a  precision  of  1.0e-6. 

Temperature  and  pressure  were  controlled  using  the  Nose-Hoover  thermostat  with  a  100  fs 
coupling  constant  and  the  Nose-Hoover  barostat  with  1000  a  fs  coupling  constant,  respectively. 
For  the  NPT  simulations,  the  three  diagonal  components  were  coupled  together  when  the 
pressure  was  computed,  and  the  dimensions  dilated  and  contracted  in  concert.  In  addition,  a  time 
step  of  1  fs  was  used  for  all  simulations  unless  otherwise  noted. 

All  simulations  were  initially  relaxed  and  annealed.  Relaxing  the  structure  involved  energy 
minimization  followed  by  four  stages  of  NVE  (constant,  particle,  volume,  and  energy) 
simulations  with  velocity  rescaling  as  the  thermostat.  The  energy  minimization  was  performed 
using  the  Polak-Ribiere  conjugate  gradient  method  with  a  maximum  of  10,000  steps.  The  energy 
and  force  tolerance  was  set  to  0  and  1.0e-8,  respectively,  with  a  maximum  of  10,000,000 
evaluations.  During  the  NVE  simulations,  the  time  step  was  increased  for  each  stage:  0.001  to 
0.010  to  0.100  to  1.000,  where  a  total  of  10,000  time  steps  was  performed  for  each  stage.  The 
velocity-rescaling  frequency,  window,  and  fraction  were  set  to  one  step,  0.01  K  and  1.0, 
respectively.  After  relaxing  the  structure,  annealing  was  performed  through  five  cycles,  where 
each  cycle  was  composed  of  50  ps  of  NPT  at  298.15  K  — »NPT  heating  at  a  rate  of  0.67  ps/K  to 
600  K— >50  ps  of  NPT  at  600  K— »NPT  cooling  at  a  rate  of  0.1656  ps/K  to  298.15  K— ^energy 
minimization.  The  energy  minimization  utilized  the  same  parameters  as  in  the  structure 
relaxation. 
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We  studied  the  effect  of  varying  equilibration  time,  system  size,  and  chain  size.  To  differentiate 
between  the  different  systems,  we  labeled  each  one  with  a  descriptor  that  combined  the  number 
of  chains  in  the  simulation  box,  number  of  monomers  per  chain,  NPT  equilibration  time  for 
reference  state,  heating  rate  of  NVT  temperature  scan,  and  additional  NVT  equilibration  time  for 
the  points  that  bracket  the  Hugoniot  point.  For  instance,  34  c-100  m-10  ns-0.015  ps/K-0, 
indicated  a  system  with  34  chains  composed  of  100  monomers,  where  10  ns  of  NPT  was  used  to 
obtain  the  reference  state.  The  NVT  temperature  scan  was  performed  using  a  heating  rate  of 
0.015  ps/K,  and  we  did  not  perform  any  additional  NVT  equilibration  time  for  the  points  that 
bracket  the  Hugoniot. 


4.  Results 


Initially,  we  considered  a  system  composed  of  34  chains  of  atactic  PMMA,  which  are  composed 
of  100  repeat  units  each  (equivalent  to  one  entanglement  per  chain).  This  equates  to 
approximately  50,000  atoms.  After  relaxing  and  annealing  the  system,  NPT  was  performed  for 
over  10  ns,  where  the  relative  values  for  the  total  energy  per  monomer  (n),  density,  temperature, 
and  pressure  vs.  time  are  shown  in  figure  2.  Relative  values  were  calculated  by  subtracting  the 
value  of  the  end  of  the  simulation  time  from  each  data  point.  From  the  figure,  we  clearly 
observed  that  the  total  energy/density  continued  to  decrease/increase  with  time  beyond  10  ns. 
The  difference  was  quite  small  O[0.1]  for  the  total  energy  per  n  and  O[0.01]  for  density),  and  the 
values  appeared  to  be  leveling  off.  It  would  be  interesting,  however,  to  note  how  equilibration 
time  could  affect  the  shape  of  the  Hugoniot  curve. 

The  Hugoniot  curve  for  34  c-100  m-10  ns-0.015  ps/K-0  is  shown  in  figure  3,  where  our 
simulated  curve  is  compared  with  experimental  data.  The  MD  results  agreed  well  with 
experiment  for  low-to-moderate  pressures  (up  to  -20  GPa),  where  deviations  were  observed  at 
higher  pressures. 

To  determine  if  a  smaller  system  would  be  able  to  reproduce  the  Hugoniot  curve,  we  also 
considered  a  system  composed  of  eight  chains  of  atactic  PMMA,  where  each  chain  was  made  up 
of  45  repeat  units  (-5000  atoms).  Again,  we  observed  that  the  total  energy  per  n  and  density 
continued  to  vary  with  simulation  time,  even  after  10  ns  of  equilibration  during  the  NPT 
simulation  under  ambient  conditions  (reference  state)  (see  figure  4). 
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Figure  2.  Relative  values  of  the  total  energy  per  monomer,  density,  temperature,  and 

pressure  vs.  time  obtained  from  NPT  simulation  for  the  reference  state  (-298  K, 
1  atm).  Relative  values  were  calculated  by  subtracting  the  value  from  the  end 
of  the  simulation  from  each  data  point.  Each  data  point  is  the  average  of  the 
last  100  ps. 
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Figure  4.  Relative  values  of  the  total  energy  per  monomer,  density,  temperature,  and 
pressure  vs.  time  obtained  from  NPT  simulation  for  the  reference  state 
(~298  K,  1  atm).  Relative  values  were  calculated  by  subtracting  the  value  from 
the  end  of  the  simulation  from  each  data  point.  Each  data  point  is  the  average 
of  the  last  100  ps. 


A  comparison  of  the  Hugoniot  curve  for  8  c-45  m-10  ns-0.015  ps/K-0  to  34  c-100  m-10  ns- 
0.015  ps/K-0  is  shown  in  figure  5.  We  observed  no  significant  difference  between  the  two 
simulated  curves,  indicating  that  systems  composed  of  shorter  and  fewer  chains  (and  fewer 
atoms)  were  adequate  for  calculating  the  Hugoniot  curve  for  PMMA.  Using  this  smaller  system, 
we  also  performed  a  study  on  the  effect  of  NPT  equilibration  time  for  the  reference  state  on  the 
Hugoniot  curve.  Instead  of  using  values  after  10  ns  of  NPT  equilibration,  we  considered  values 
at  2  ns  and  100  ps.  Hugoniot  curves  for  these  two  systems  are  shown  in  figures  6  and  7,  where 
8  c-45  m-10  ns-0.015  ps/K-0  did  not  noticeably  vary  from  8  c-45  m-2  ns-0.015  ps/K-0  and 
8  c-45  m-100  ps-0.005  ps/K-0.  From  this,  we  inferred  that  the  calculation  was  relatively 
insensitive  to  the  small  change  in  the  reference  total  energy,  density,  pressure,  and  temperature 
achieved  by  performing  the  NPT  simulation. 
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Figure  5.  Hugoniot  curve  for  8  c-45  m-10  ns-0.015  ps/K-0  and  34  c-100  m-10  ns-0.015  ps/K-0. 


Figure  6.  Hugoniot  curve  for  8  c-45  m-2  ns-0.015  fis/K-0  and  8  c-45  m-10  ns-0.015  fis/K-0. 


Figure  7.  Hugoniot  curve  for  8  c-45  m-100  ps-0.005  gs/K-0  and  8  c-45  m-2  ns-0.015  ps/K-0. 


To  determine  whether  the  method  was  sensitive  to  the  number  of  monomers  per  chain,  we 
considered  chains  composed  of  either  one  (360  c-1  m-100  ps-0.005  ps/K-0,  -5000  atoms)  or  four 
(33  c-4  m-100  ps-0.005  ps/K-0,  -2000  atoms)  monomers.  The  Hugoniot  curves  for  these  two 
simulations  are  shown  in  figures  8  and  9.  A  clear  difference  between  360  c-1  m-100  ps-0.005 
ps/K-0  and  8  c-45  m-100  ps-0.005  ps/K-0  was  shown,  indicating  that  the  concept  of  chain 
connectivity  was  important  for  the  Hugoniot  calculation.  Increasing  the  chain  to  four  monomers, 
the  calculated  Hugoniot  curve  agreed  well  with  the  curve  that  was  produced  with  a  chain 
composed  of  45  monomers. 


Particle  Velocity  (km/s)  V/Vo 


Figure  8.  Flugoniot  curve  for  8  c-45  m-100  ps-0.005  ps/K-0  and  360  c-1  m-100  ps-0.005  ps/K-0. 


Particle  Velocity  (km/s)  V/Vo 


Figure  9.  Hugoniot  curve  for  8  c-45  m-100  ps-0.005  ps/K-0  and  33  c-4  m-100  ps-0.005  ps/K-0. 
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This  observation  allowed  us  to  look  towards  quantum  mechanics  to  calculate  Hugoniot  curves. 
Quantum  can  be  used  to  look  at  chain  breaking  which  occurs  at  high  pressures,  but  this  is 
computationally  expensive.  As  a  result,  only  small  systems  where  the  polymers  can  be 
represented  by  a  few  monomers  are  computationally  feasible. 

To  determine  whether  it  was  necessary  to  further  equilibrate  the  points  around  the  Hugoniot, 
which  were  used  to  interpolate  the  Hugoniot  temperature  and  pressure,  we  performed  an  additional 
10  ns  of  NVT  after  the  temperature  scan.  Performing  this  additional  equilibration  did  decrease 
the  total  energy  and  pressure,  as  seen  in  figures  10  and  1 1,  for  34  c-100  m-10  ns-0.015  ps/K-10  ns 
and  8  c-45  m-10  ns-0.015  ps/K-10  ns,  respectively.  However,  it  did  not  decrease  significantly 
where  values  changed  less  than  0(1),  even  after  an  additional  10  ns  of  equilibration. 
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Figure  10.  Relative  values  of  the  total  energy  per  monomer,  density,  temperature, 

and  pressure  vs.  time  obtained  from  NVT  simulation  at  elevated  pressure 
and  temperature  (-1280  K,  28  GPa).  Relative  values  were  calculated  by 
subtracting  the  value  from  the  end  of  the  simulation  from  each  data  point. 
Each  data  point  is  the  average  of  the  last  100  ps. 
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Figure  1 1 .  Relative  values  of  the  total  energy  per  monomer,  density,  temperature, 
and  pressure  vs.  time  obtained  from  NVT  simulation  at  elevated 
pressure  and  temperature  (-1300  K,  28  GPa).  Relative  values  were 
calculated  by  subtracting  the  value  from  the  end  of  the  simulation  from 
each  data  point.  Each  data  point  is  the  average  of  the  last  100  ps. 


The  Hugoniot  curves  calculated  with  and  without  this  additional  10-ns  NVT  equilibration  are 
shown  in  figures  12  and  13  for  34  c-100  m-10  ns-0.015  ps/K  and  8  c-45m-10  ns-0.015  ps/K. 
There  were  no  significant  differences  between  the  two  curves,  thus  implying  that  this  additional 
equilibration  was  not  necessary  for  PMMA. 


Figure  12.  Hugoniot  curve  for  34  c-100  m-10  ns-0.015  ps/K-10  ns  and  34  c-100  m-10  ns-0.015  ps/K-0. 
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Figure  13.  Hugoniot  curve  for  8  c-45  m-10  ns-0.015  gs/K-10  ns  and  8  c-45  m-10  ns-0.015  gs/K-0. 


5.  Conclusion 


Through  MD  simulations  implemented  through  the  open  source  code  LAMMPS,  PCFF  was 
accessed  for  predicting  PMMA  subjected  to  extreme  pressures  and  temperatures.  In  addition,  we 
also  studied  the  effect  of  simulation  parameters  such  as  equilibration  time,  system  size,  and  chain 
size  on  the  Hugoniot  curve.  For  the  former,  we  found  good  agreement  between  the  experimental 
and  simulated  Hugoniot  curves  at  low  and  moderate  pressures  and  temperatures.  For  the  latter, 
we  found  that  the  NPT  equilibration  time  for  the  reference  state  had  a  minimal  effect  on  the 
Hugoniot  curve  since  annealing  and  relaxing  the  system  adequately  optimized  the  system.  In 
addition,  we  found  that  the  Hugoniot  points  could  be  determined  by  a  quick  temperature  scan, 
where  further  equilibration  of  the  points  that  bracket  the  Hugoniot  point  was  unnecessary. 
Finally,  we  also  found  the  smaller  system  sizes  were  adequate  for  Hugoniot  simulations. 
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